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ABSTRACT 


Goals  of  this  research  have  been  to  identify  physical  processes  that  determine  the  dynamics 
of  the  marine  cloud  layers  and  to  quantify  roles  of  turbulence,  convection  and  thermal  radiation  that 
play  in  formation,  dissipation  and  stability  of  the  marine  cloud  layers.  And  immediate  objectives  of 
the  research  are  to  advance  turbulence  models,  use  efficient  numerical  schemes,  develop  computer 
simulation  programs,  simulate  the  marine  cloud  layers  and  compare  computer  results  with  published 
experimental  data  on  the  marine  cloud  layers  so  as  to  yield  insights  into  the  cloud’s  physical 
processes. 

For  these  objectives,  two  theoretical  models,  using  the  second-order-turbulence  closure  and 
the  large-eddy  simulation  (LES),  respectively,  have  been  developed.  In  addition,  efforts  have  been 
made  to  develop  a  hybrid  model  that  is  based  upon  a  framework  of  the  LES  model  but  uses 
turbulence  values  predicted  by  a  second-order-closure  model.  This  hybrid  model  preserves  details 
of  turbulence  but  eliminates  the  need  of  continuous  evaluation  of  the  multi  dimensional  turbulence 
equations.  The  model  will  improve  simulation  efficiency  and  conserve  computational  resources.  An 
evaluation  of  the  model  has  shown  excellent  agreement  between  simulation  results  with  relevant 
experimental  data  retrieved  from  reliable  web  site  resources. 

In  this  report,  a  progressive  development  of  the  second-order-closure,  LES  and  hybrid 
turbulence  models  for  simulation  of  the  marine  cloud  layers  is  described  and  a  comparison  of 
theoretical  results  with  experimental  data  is  presented  to  yield  better  insights  into  the  cloud  s  physical 


processes. 
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1.  INTRODUCTION 


Goals  of  this  research  have  been  to  identify  physical  processes  that  determine  the  dynamics 
of  the  marine  cloud  layers  and  to  quantify  roles  of  turbulence,  convection  and  thermal  radiation  that 
play  in  formation,  dissipation  and  stability  of  the  marine  cloud  layers.  And  immediate  objectives  of 
the  research  are  to  advance  turbulence  models,  use  efficient  numerical  schemes,  develop  computer 
simulation  programs,  simulate  the  marine  cloud  layers  and  compare  theoretical  results  with  published 
experimental  data  on  the  marine  cloud  layers  so  as  to  yield  insights  into  the  cloud’s  physical 
processes. 

Studying  the  multi  dimensional  cloud  layers  requires  a  turbulence  model  that  permits  high- 
resolution  simulation  of  turbulence  of  different  scales  and  domain  sizes.  For  high  resolution,  detailed 
second-order-closure  model  of  turbulence  may  be  used,  and  for  calculation  efficiency,  the  large-eddy- 
simulation  (LES)  model  may  be  developed.  Two  computer  simulation  programs,  using  the  second- 
order-closure  model  (Chi,  1996  and  1998a)  and  the  LES  model  (Chi,  1998b),  respectively,  have  been 
developed  under  this  research  effort.  Additional  efforts  have  been  made  to  develop  a  hybrid  model 
that  is  based  upon  a  framework  of  the  LES  model  but  uses  turbulence  values  predicted  by  a  second- 
order-closure  model.  In  addition,  confidence  in  the  simulation  model  is  established  by  comparing 
simulation  results  with  experimental  data  retrieved  from  reliable  web  site  resources. 

In  this  report,  a  progressive  development  of  the  second-order-closure,  LES  and  hybrid 
turbulence  models  for  simulation  of  the  marine  cloud  layers  will  be  documented  first.  Numerical 
experiments  will  be  carried  out  to  examine  calculation  efficiencies  and  characteristics  of  the  marine 
cloud  layers.  Experience  in  tracking  relevant  experimental  data  from  reliable  web  site  resources  will 
be  described.  To  increase  insight  into  the  marine  cloud  layers,  the  computer  simulation  results  are 


compared  with  experimental  data. 

2.  THEORETICAL  ■  MODELS  FOR  THE  MARINE  £LQ1 JR  LAYERS 

For  high  resolution,  a  detailed  second-order-closure  model  of  turbulence  has  been 
developed.  For  calculation  efficiency,  a  LES  model  has  been  developed.  Details  of  these  models 
have  been  reported  in  several  papers  by  the  author  (Chi,  1996,  1998a,  and  1998b).  For  ready 
references,  reprints  of  these  three  papers  are  appended  to  the  report.  In  addition,  a  hybrid 
turbulence  model  that  is  built  upon  a  framework  of  the  LES  model  but  uses  turbulence  values 
predicted  by  a  second-order-closure  model  has  been  used  to  predict  stability  of  multi  dimensional 
marine  cloud  layers.  This  hybrid  model  preserves  details  of  turbulence  but  eliminates  the  need  of 
evaluating  the  multi  dimensional  turbulence  equations,  and  It  improves  prediction  efficiency  and 
conserve  computational  resources.  In  addition,  prediction  results  have  been  shown  in  excellent 
agreement  with  relevant  experimental  data  retrieved  from  reliable  web  site  resources.  Presented 
below  in  Section  2.1  is  a  second-order-closure  model  for  calculating  the  turbulence  values  and  in 
Section  2.2  is  the  framework  of  a  LES  model.  In  Section  2.3,  a  hybrid  model  that  is  based  upon  the 
LES  framework  but  uses  the  second-order-closure  turbulence  values  is  described  to  predict  stability 
of  the  marine  cloud  layers. 

2.1  A  Second-Order-Closure  Model  of  Turbulence 

When  these  assumptions  are  made:  (1)  vertical  coordinate  z  is  in  the  dominant  direction  of 
turbulent  mixing  of  atmospheric  air,  (2)  at  far  above  a  sea  surface,  the  geostrophic  balance  is 
maintained,  (3)  velocity  and  temperature  of  vapor  and  liquid  moisture  are  in  equilibrium,  and  (4) 
Boussinesq  approximations  are  used,  the  conservation  equations  for  momentum,  enthalpy  and  total 
moisture  of  atmospheric  air  can  be  written  as: 
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where  the  upper-case  dependent  variables  represent  the  mean  fields  and  the  corresponding  lower¬ 
case  variables  represent  turbulent  fluctuations  of  the  same  variables.  Thus,  U,  V  and  W  are  mean 
velocity  components  in  the  x,  y  and  z  directions  and  at  the  time  t.  (Ug,Vg)  are  geostrophic  wind 
components  in  (x,y)  directions.  0  is  the  mean  moist-air  potential  temperature  which  is  defined 
as  [T-T0+  (gz  +  LQ^/Cp],  and  Q  the  mean  total  moisture  mixing  ratio  which  is  defined  as  (Qv+Q,). 
L  is  the  water  latent  heat  of  vaporization,  Cp  the  constant  pressure  specific  heat,  T  the 
temperature,  T0  the  referenced  temperature  at  z  equal  to  zero,  the  water  vapor  mixing  ratio, 
Q,  the  liquid  water  mixing  ratio,  g  the  gravitational  acceleration,  f  the  Coriolis  parameter,  ®v  the 
virtual  dry  potential  temperature  which  is  defined  as  [T(l  +  1.6O90v-Q)-To+gz/Cp],  g  the  standard 
density,  v  the  kinematic  viscosity,  and  P  the  buoyancy  coefficient,  The  second-order  correlations 
in  these  equations  represent  the  mean  turbulent  fluxes  of  momentum,  enthalpy  and  moisture 
fluxes. 

Using  the  second-order-closure  assumption  (Mellor  and  Yamada,  1974  and  Moeng  and 
Arakawa,  1980),  transport  equations  for  calculating  the  turbulent  fluxes  may  be  written  as  follows: 
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In  equations  5  through  18,  the  time-derivative  terms  on  the  left-hand  side  model  the  transient 
variation  of  turbulence  correlations.  The  second-order  derivative  terms  on  the  right-hand  side  model 
turbulent  diftusion.  While  the  production  of  turbulence  due  to  buoyancy  is  modeled  by  terms  with 
buoyant  coefficient  P,  the  production  of  turbulence  due  to  friction  is  modeled  by  products  of  second- 
order  correlations  and  gradients  of  mean  variables.  Terms  with  coefficients  B,  C,  E  and  F  represent 
the  turbulent  redistribution.  The  characteristic  length  scale  X  is  equal  to  the  value  of  the  Blackadar's 
or  the  diffusion-length  scale  -  XB  or  XD  -  whichever  is  the  smallest.  Turbulent  dissipation  is  modeled 
by  terms  with  coefficient  D.  Coefficients  k.  A,  B,  C,  D,  E  and  F  have  been  determined  semi- 
empirically  (Chi,  1994),  having  the  values  equal  to  0.35,  0.21,  0.46,  0.053,  0.132,  0.44,  and  0.23, 
respectively.  Numerical  procedures  and  computer  programs  have  been  developed  to  solve 
equations  1  through  18  with  appropriate  boundary  conditions  (Chi,  1996  and  1998a).  Prediction 
results  using  the  second-order-closure  model  will  be  presented  and  discussed  in  Section  4. 1  of  this 
report. 
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2.2  A  LES  Model  for  Marine  Cloud  Layers 

When  these  assumptions  are  made:  (1)  the  Coriolis  force  is  negligible,  (2)  velocity  and 
temperature  of  vapor  and  liquid  moisture  are  in  equilibrium,  and  (3)  Boussinesq  approximations 
are  used,  the  conservation  equations  for  momentum,  enthalpy,  total  moisture  and  turbulent  energy 


of  atmospheric  air  can  be  written  as: 
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In  above  equations,  E  is  the  turbulent  energy  defined  as  q2/2,  and  the  buoyancy  terms  associated  with 
the  virtual  dry  potential  temperature  0V  have  been  defined  as  follows: 
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The  rate  of  dissipation  within  the  grid  volume  and  the  subgrid  eddy  coefficient  may  be  parameterized 

• 

through 
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and  the  eddy  diffusivity  coefficients  have  been  defined  as  follows: 
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And  the  turbulent  Prandti  numbers  for  enthalpy  and  moisture  diflusion  (o„  and  oj  are  both  equal  to 
0.75.  Numerical  procedures  and  computer  programs  have  been  developed  to  solve  equations  19 
tough  31  with  appropriate  boundary  conditions  (Chi,  1998b).  Prediction  results  using  the  LES 

model  will  be  presented  and  discussed  Section  4.2  of  this  report. 

In  Section  2.3,  a  hybrid  model  that  retains  high  resolution  with  the  second-order-closure 
model  described  in  Section  2.1  and  calculation  efficiency  with  the  LES  model  described  in  Section 


2.2  will  now  be  presented. 

2.3  A  Hybrid  Turbulence  M»Hp1  for  Marine  Cloud  Layers 

Second-  and  higher-order  turbulent  models  have  succeeded  in  advancing  theoretical 


understanding  of  the  marine  cloud  layers.  Complexity  of  those  models  often  makes  long-term 
simulation  of  the  cloud  layer  over  an  extensive  period  prohibitively  expensive.  A  hybrid  model  that 
retains  calculation  efficiency  of  the  LES  model  and  turbulence  details  of  the  second-order-closure 


model  has  been  developed.  The  hybrid  model  uses  Equations  19  through  23  from  the  LES  model  for 
conservations  of  momentum,  enthalpy,  moisture  and  mass,  respectively: 
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•  But,  instead  of  the  turbulent-kinetic-energy  Equation  24,  the  second-order-closure  model  equations 

1  through  18  will  be  used  in  evaluating  the  eddy-diffusivity-coefficient  K  values  defined  by  Equations 
29  through  32.  Numerical  Procedure  and  computer  algorithms  (described  in  Section  3  below)  have 
^  been  developed  to  evaluate  the  hybrid  model.  Prediction  results  using  the  hybrid  model  will  be 

presented  in  Section  4.3  and  verified  by  comparison  with  experimental  data  in  Section  5.2. 

0  3.  NUMERICAL  procedure  and  computer  algorithm 

The  hybrid  model  presented  in  Section  2.3  uses  the  turbulence  equations  1  through  18,  the 
conservation  equations  19  through  23,  and  the  diffusivity  coefficients  defined  by  Equations  29 
^  through  32.  Numerical  procedures  and  computer  algorithms  have  been  developed  to  simulate  the 

marine  cloud  layers.  It  can  be  observed  above  that  conservation  and  turbulence  transport  equations 
f  1  through  23  described  above  can  be  written  in  the  following  general  form: 
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The  boundary  and  initial  conditions  may  be  written  as: 
9a 
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a(x,z,t=0)=aT](x,z),  At  initial  time  zero 


(39) 


In  above  equations,  a(x,z,t)  values  are  variable  values  of  velocity  components  U  and  W,  potential 
temperature  0,  and  moisture  mixing  ratio  Q;  <|)(x,z,t)  values  are  the  source/sink  strengths  for  c c. 
Items  in  equation  37  describe  the  variables’  time  derivative,  convection,  diffusion  and  source/sink 
strength,  respectively,  'a'  and  'b'  values  in  equation  38  are  used  to  define  the  appropriate  boundary 
conditions. 

A  finite-volume  difference  scheme  (Patanka,  1980)  was  used  for  descritization  of  variables. 
The  flow  domain  to  be  simulated  was  divided  into  small  rectangular  elements;  shown  in  a  figure 
below  are  examples  of  several  such  elements.  It  can  be  seen  in  the  figure  that  variable  values  at  node 
points  of  the  elements  are  strategically  located.  The  node  point  for  velocity  value  is  located  in  the 
middle  of  the  rectangular  edge  and  the  node  point  for  other  variables  is  located  in  the  center  of  the 
rectangular  box. 

A  Figure  Showing  an  Element  with  Neighboring  Variables 
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Using  the  velocity  values  and  values  of  a  shown  in  the  figure  and  employing  an  upwind  logic, 
the  conservation  law  may  be  applied  to  obtain  an  expression  for  the  value  of  a  at  the  node  point  P 
in  terms  of  the  a-values  at  its  neighboring  node  points: 


apap=aEaE+an^w+a7aT+aBaB+b 

(40) 

where  coefficient  values  ‘a’  and  ‘b’  may  be  expressed  in  terms  of  the  diffusion  and  flux  parameters 

defined  as  follows: 
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In  addition,  in  solving  U  and  W  values  using  equation  40,  the  initial  pressure  values  will  have  to 
be  the  estimated  values;  consequently,  the  equation  40  will  yield  initially  approximate  Uw  and  Ww 
values.  Improved  U  and  W  values  may  be  calculated  from  their  approximate  values  using  equations: 
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Substituting  U  and  W  values  in  equations  47  to  50  into  the  continuity  equation  23  yields  a  set  of 
linear  equations  for  pressure  values  at  the  solution  domain  s  node  points,  they  may  be  used  to  solve 
for  improved  pressure  values.  So  the  process  may  be  repeated  to  iterate  alternatively  for  improved 

values  of  velocity,  pressure  and  other  variables  at  the  node  points. 

Using  the  linear  equations  discussed  above  for  variable  values  U,  W,  0,  Q, E  and  P  at  node 
points,  a  computer  program  has  been  written  to  simulate  dynamics  of  the  marine  cloud  layers. 

4.  COMPUTER  STMTTT  ATTON  RESULTS 

4.1  Sprnnri-Order-Closure  Model  Simulation  Results 

The  second-order-closure  model  of  turbulence  described  in  Section  2.1  has  been  used  to 
predict  transient  exchange  of  moisture  in  the  sea/air  interface,  mixing  of  moisture  in  the  marine 
atmosphere,  and  formation  of  the  cloud  layer.  Graphs  plotted  in  figure  1  were  reported  in  a  previous 
paper  written  by  the  author  (Chi,  1996);  they  were  used  as  initial  conditions  for  this  study. 

The  initial  conditions  shown  in  figure  1  were  obtained  on  the  assumptions  that  the  geostrophic 
wind  in  the  x-direction  was  at  lOm/s,  and  Coriolis  parameter  f  is  equal  to  l-OxlO"4  s’1,  the  potential 
temperature  was  chosen  to  be  well-mixed  at  290K,  and  the  total  moisture  mixing  ratio  corresponding 
to  100  percent  relative  humidity  at  top  of  the  marine  planetary  boundary  layer  (MPBL).  Just  above 
the  MPBL  top  at  height  equal  to  1  kilometer,  there  was  an  inversion  layer  of  400-m  thick,  in  which 


# 


# 


* 
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Fig.  1 :  Initial  Velocity-Components  and  Potential-Temperature  Values 

temperature  was  increasing  at  0.2  K/m  and  the  moisture  mixing  ratio  was  decreasing  at  0.008 
g/kg.m,  respectively.  The  sea-surface  temperature  was  allowed  to  vary  diumally  within  the  range  of 

# 

280  to  290F.  Conditions  shown  in  figure  2  were  for  the  instance  when  the  sea  surface  was  at  290  K. 

%  For  the  simulation  run  using  initial  conditions  shown  in  figure  1,  it  was  assumed  that  the  sea 

# 
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surface  temperature  was  raised  abruptly  to  and  then  maintained  at  298  K,  and  the  air  total  moisture 
ratio  at  the  sea  surface  was  maintained  at  the  saturation  state.  After  simulation  runs  over  a  period 
of  forty  hours,  a  large  body  of  data  was  generated.  Plotted  in  figures  2  are  predicted  contours  of 
temporal  mean  physical-property  values  of  air  in  the  marine  atmosphere;  plotted  in  figures  3  are 
predicted  contours  of  Reynolds-stress  and  turbulent-flux  values  in  the  marine  atmosphere.  Many 
interesting  characteristics  of  mixing  in  the  simulation  domain  can  be  observed  in  these  contours. 
Firstly,  rapid  interaction  at  the  air/sea  interface  can  be  observed  during  the  initial  period  of  zero  to 
600  minutes.  It  is  followed  by  a  calmer  development  of  the  marine  planetary  boundary  layer  for  about 
600  minutes.  Then,  during  the  next  600  minutes  transfer  of  enthalpy  and  moisture  continues,  as  can 
be  observed  from  the  predicted  contours  of  the  turbulent  thermal  and  moisture-flux  values  shown  in 
figures  2C  and  2D.  Also,  can  be  observed  in  Fig.  2D  is  the  formation  of  cloud  starting  at  the  200th 
minute,  rapid  deepening  of  the  cloud  layer  during  the  second  interval  of  200  to  1200  minutes,  and 
slower  growth  of  the  cloud  layer  during  the  period  of  1200  to  1800  minutes.  Finally,  a  steady  state 
is  established  at  around  the  2400th  minute. 

Plotted  in  figures  4  and  5  are  the  predicted  steady-state  conditions  at  the  end  of  the  simulation 
period  of  the  40th  hour.  The  graphs  shown  in  figure  4  can  be  compared  with  those  in  figure  1  for  the 
initial  conditions.  Changes  that  have  been  made  during  these  forty  hours  can  be  observed.  Warm  and 
moist  air  at  a  sea  surface  has  resulted  in  an  unstable  mixing  layer.  Variations  of  both  the  mean- 
quantity  and  the  turbulent-flux  profiles  can  be  observed  in  figure  4.  Specifically,  if  mean  horizontal 
velocity  components  plotted  in  figures  1 A  and  turbulent  kinetic  energy  values  plotted  in  figure  IB 
are  compared  with  those  corresponding  values  plotted  in  figure  4A  and  4H,  thickening  of  the 
unstable  boundary  layer  due  to  intense  turbulent  exchange  can  be  observed  in  figure  4A  and  4H. 


Fig.  2:  Predicted  Contours  of  Mean  Atmospheric  Quantities 
2A  -  Mean  Wind  Velocity,  m/s,  2B  -  Mean  Potential  Temperature,  K, 

2C  -  Mean  total  moisture  mixing  ratio,  g/kg,  2D  -  Mean  liquid  moisture  mixing  ratio,  g/kg 


Fig.  3:  Predicted  Contours  of  Turbulent  Fluxes 
3A  -  Turbulent  kinetic  energy,  mA2/sA2,  3B  -  Principal  Reynolds  stress,  mA2/sA2 
3C  -  Turbulent  thermal  flax,  m.K/s,  3D  -  Turbulent  moisture  flux,  m.kg/kg.sA3 
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Fig.  4:  Predicted  Mean  Quantities  and  Turbulent  Fluxes 

4A  -  Mean  velocity;  4B  -  Reynolds  stress;  Q 

4C  -  Potential  temperature;  4D  -  Thermal  flux; 

4E  -  Total  moisture  ratio;  4F  -  Moisture  flux; 

4G  -  Liquid  moisture  ratio;  4H  -  Turbulent  energy. 
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Fig.  5:  Initial  steady  state  at  the  zeroth  hour 


Fig.  6:  Contours  of  superposed  stream  functions  with 
entrainment  at  the  top 


4.2  LES  Model  Simulation  Results 


The  LES  model  presented  in  Section  2.2  has  been  used  to  study  stability  of  the  marine  cloud 
layers.  Graphs  plotted  in  figure  4  and  5  predicted  by  a  second-order-closure  model  presented  in 
Section  4. 1  above  can  be  used  as  a  set  of  initial  values  in  this  numerical  experiment.  It  may  be  noted 
that  solutions  for  this  set  of  graphs  were  obtained  by  using  the  main  stream  wind  velocity  U  equal  to 
10  m/s.  To  investigate  effects  of  cloud-top-warm-air  entrainment  on  stability  of  the  cloud  layer 
shown  in  figure  5,  a  stream  function  shown  in  figure  6  will  be  superimposed  to  the  main  flow.  In 
generating  the  stream  function  shown  in  figure  6,  the  top-down  entrainment  was  assumed  to  be  at  a 
maximum  rate  of  6  cm/s  is  at  the  top-left  comer  and  the  rate  was  reduced  sinusoidally  to  zero  at  the 
top-right  comer.  In  addition,  it  is  assumed  that  potential  temperature  of  the  entrained  air  is  at  five 

degrees  centigrade  higher  than  that  of  the  cloud  at  the  top. 

Using  the  initial  conditions  and  entrainment  rates  described  above,  a  LED  simulation  computer 
program  has  been  run.  Shown  in  figure  5  are  contours  of  the  initial  steady-state  potential  temperature 
values,  total  moisture  mixing  ratio  values  and  liquid  moisture  Mixing  values,  respectively.  Predicted 
dynamic  responses  of  the  cloud  layer's  liquid  moisture  content  to  the  warm-air  entrainment  at  the 
top  are  shown  in  figure  7.  From  snapshots  shown  in  this  figure  of  the  cloud  contours  at  different 
times  (i.e,  at  half,  one,  five  and  ten  hours  from  the  start  of  the  simulation  run),  dissipation  of  the 

cloud  layer  can  be  observed. 

A  comparison  of  this  example  with  that  presented  in  Section  4. 1,  it  can  be  observed  that  the 
LES  model  is  superior  to  the  second-order-closure  model  in  calculation  efficiency.  However,  details 
of  turbulence  are  lost  in  trading  for  this  efficiency. 
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4.3  Hybrid  Mo  dpi  Simulation  Results 

In  order  to  preserve  details  of  the  second-order-closure  turbulence  and  calculation  efficiencies 
of  large-eddy  simulation,  a  computer  program  using  the  hybrid  model  developed  in  section  2.3  has 
been  written  to  study  stability  of  the  marine  cloud  layer.  Numerous  runs  of  the  computer  program 
have  been  made  to  simulate  dynamic  responses  of  marine  cloud  layers  under  a  variety  of  conditions. 
To  facilitate  later  comparison  of  simulation  results  with  TOGA-COARE  data  (Tropical  Ocean  Global 
Atmosphere  -  Coupled  Ocean  Atmosphere  Response  Experiment)  documented  by  Tao  and  Simpson 
(1993),  Miller  and  Riddle  (1994)  and  Lin  and  Johnson  (1996),  computer  runs  have  been  made  with 
internet  retrieved  experimental  stream-function  and  boundary-condition  values  as  input  data.  As  an 
example,  figure  8  shows  an  initial  experimental  dataset  for  the  stream-function  values,  and  figure  9 
shows  snapshots  of  cloud  profiles  predicted  by  the  hybrid  model,  at  several  different  instances, 
starting  from  the  zeroth  hour  to  the  24th  hour. 


Fig.  8:  Contours  of  a  set  of  retrieved  experi  rental  stream  function  values 
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Fig.  9:  Snapshots  of  cloud  contours  predicted  by  a  hybrid  model  of  turbulence 
(Starting  from  an  initial  experimental  dataset  at  the  zeroth  hour) 


5.  DATA  TRACKING  AND  THEORETICAL  VERIFICATION 

5.1  Internet  Atmospheric  Data  Trackin2 

# 

Several  data  sources  such  as  Goddard  Cumulus  Ensemble  (GCE),  GDACC  (Goddard 
Distributive  Active  Archive  Center),  FIFI  (First  ISLCP  Field  Experiment)  and  TOGA-COARE 
(Tropical  Ocean  Global  Atmosphere  -  Coupled  Atmosphere  Response  Experiment)  were  reviewed  % 

to  determine  the  time  interval,  grid  size  and  altitude  criteria.  It  was  determined  that  we  should  secure 
data  with  time  and  grid  intervals  at  six  hours  and  one  mile,  respectively.  Two  sources  of  data  TOGA- 
COARE  (Webster  and  Lukas,  1992)  and  GCE  (Tao,  1993)  were  then  considered. 

TOGA-COARF  DATA:  The  scientific  goals  of  COARE  are  to  describe  and  understand:  (1) 
the  principal  processes  responsible  for  the  coupling  of  the  ocean  and  the  atmosphere  in  the  western  0 

Pacific  warm-pool  systems,  (2)  the  principal  atmospheric  processes  that  organize  convection  in  the 
warm-pool  region,  (3)  the  oceanic  response  to  combined  buoyancy  and  wind-stress  forcing  in  the 

.  m 

western  Pacific  warm-  pool  region,  and  (4)  the  multiple-scale  interactions  that  extend  the  oceanic  and 
atmospheric  influence  of  the  western  Pacific  warm-pool  system  to  other  regions  and  vice  versa.  To 
carry  out  the  goals  of  TOGA  COARE,  three  components  of  a  major  field  experiment  have  been  ^ 

defined:  interface,  atmospheric,  and  oceanographic.  The  experimental  design  calls  for  a  complex  set 
of  oceanographic  and  meteorological  observations  from  a  variety  of  platforms  that  carry  out  remote 
and  in  situ  measurements.  The  resulting  high-quality  dataset  is  required  for  the  calculation  of  the 
interfacial  fluxes  of  heat,  momentum  and  moisture,  and  to  provide  ground  thruth  for  a  wide  range  of 
remotely  sensed  variables  for  the  calibration  of  satellite-derived  algorithms.  The  ultimate  objective  ^ 

of  the  COARE  dataset  is  to  improve  air-sea  interaction  and  boundary-layer  parameterizations  in 
models  of  the  ocean  and  the  atmosphere,  and  to  validate  coupled  models. 
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Internet  web-site  data  was  used  to  review  the  data  sets  of  TOGA-COARE.  Several  sites  for 
data  were  reviewed.  The  web  site,  http://kiwi.atmos.colostate.edu/scni/toga-coare.html  was  found 
to  contain  mean  data  over  the  TOGA-COARE  IFA  region,  which  were  thought  appropriate  for 
testing  the  model.  Using  FTP  commands  the  data  files  were  transferred  to  a  main  frame  computer. 
Using  uncompress  command  the  compressed  files  were  uncompressed  and  converted  to  ASCII  files 
and  transmitted  to  a  UDC  workstation  via  electronic  mails.  However,  owing  to  their  coarse  grid  sizes, 
higher  resolution  datasets  are  required  for  the  model  refinement. 

GTE  MODELS  DATA:  The  Goddard  Cumulus  Ensemble  Model  (GCE)  is  maintained  by 
Mesoscale  Atmospheric  Processes  Branch  (MAPB)  at  Goddard  Space  Flight  Center  (NASA/GSFC). 
Scientists  in  NASA/GSFC  continuously  maintain  and  upgrade  the  system  (Tao  and  Simpson,  1993), 
and  the  GCE  model  is  simulated  and  operated  by  GSFS  system  analysts  and  operators  (Private 
communications).  The  GEC  model  may  be  used  to  interpolate  experimental  data  with  great  degrees 
of  accuracy.  The  following  procedures  were  followed  to  obtain  the  data  from  both  sources. 

Super  computer  CRAY  at  NASA  was  used  for  obtaining,  organizing,  processing  the  data. 
The  data  files  obtained  in  this  process  were  transmitted  to  the  UDC  workstation  via  electronic  mail. 
For  GCE  Model  data,  scientists  at  GSFC  provided  the  author  and  his  co-workers  at  UDC  with  the 
database  and  a  Fortran  program  including  a  list  of  parameters  via  electronic  mail.  Fortran  programs 
were  further  developed  and  modified  the  data  into  appropriate  formats  for  the  files.  The  data  from 
GCE  model  were  processed  using  the  Fortran  programs  in  UNIX  environment.  Then.  They  were 
transmitted  to  the  UDC  workstation  in  two  stages. 

The  First-stage  data  transmitted  included  the  following  parameters: 
parameter  (nx=512,nz=22) 
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parameter  (nhr=48)  !  number  of  hours  model  simulations  (h) 

parameter  (interval=2 1600)  !  Interval  to  read  in  data  (s) 

parameter  (ntimes=nhr  *3 600/interval) 


real  latent  (nx,ntimes) 
real  pO  (nz) 
real  qcg(nx,nz,ntimes) 
real  qci(nx,nz,ntimes) 
real  qcl(nx,nz,ntimes) 
real  qcs(nx,nz,ntimes) 
real  qm(nx,nz,ntimes) 
real  rad_lw(nx,nz,ntimes) 
real  rad_sw(nx,nz,ntimes) 
real  sensible(nx,ntimes) 
real  temp(nx,nz,ntimes) 
real  u(nx,nz,ntimes) 
real  v(nx,nz,ntimes) 
real  vap(nx,nz,ntimes) 
real  w(nx,nz,ntimes) 
real  zl(nz) 
real  z2(nz) 


!  surface  latent  heat  flux  (W/m**2) 

!  air  pressure  (mb) 

!  graupel  mixing  ratio  (g/kg) 

!  cloud  ice  mixing  ratio  (g/kg) 

!  cloud  water  mixing  ratio  (g/kg) 

!  snow  mixing  ratio  (g/kg) 

!  rain  water  mixing  ratio  (g/kg) 

!  long  wave  radiative  heating  rate  (K/hr) 

!  short  wave  radiative  heating  rate  (K/hr) 

!  surface  sensible  heat  flux  (W/m**2) 

!  air  temperature  (K) 

!  u-wind  speed  (m/s) 

!  v-wind  speed  (m/s) 

!  water  vapor  mixing  ratio  (g/kg) 

!  w-wind  speed  (m/'s) 

!  model  grid  height  at  zl  levels  -  for  most  variables  here  (m) 
!  model  grid  height  at  z2  levels  -  where  2  is  (m) 


nx=5 12  (512  miles,  data  grid  =  one  mile) 
nz=22  ->  pressure  levels 
ntimes=32  (6-hrs  intervals) 


The  Second-stage  data  transmitted  included  the  following  parameter: 


real  tke(nx,nz) 

real  du(nx,nz) 

real  dv(nx,nz) 

real  dw(nx,nz) 

nz=22  ->  pressure  levels 

ntimes=32  (6-hrs  intervals) 


!  turbulent  kinetic  energy  (cm*cm/s/s) 

!  u-momentum  turbulence  rate  (cm/s/s) 

!  v-momentum  turbulence  rate  (cm/s/s) 

!  w-momentum  turbulence  rate  (cm/s/s) 


The  database  obtained  from  the  GCE  model  was  large.  Eight  six-hours  data  were  organized  in  eight 


different  data  files,  each  file  containing  one  six-hour  interval  of  data.  The  data  files  were  transmitted 


via  electronic  mails. 


DATA  PRESENTATION:  For  graphic  presentation,  they  were  further  decoded  at  UDC 
by  several  Fortran  programs,  and  were  translated  into  formats  recognizable  by  data  plotting  software 
packages  -  Surfer  and  Grapher  -  which  are  available  in  the  author’s  laboratory.  Presented  below  are 
plotted  sample  data  values  over  a  domain  with  a  height  of  3 -km  and  a  horizontal  distance  of  500-km, 
and  over  a  time  period  of  24  hours  at  12-hours  intervals: 

•  Contours  of  empirical  wind  flow  stream  functions  values 

•  Contours  of  empirical  total  moisture  mixing  ratio  values 

•  Contours  of  empirical  liquid  moisture  mixing  ratio  values 

•  Sea/air  interface  sensible  and  latent  heat  flux  values 


Fig.  10:  Contours  of  experimental  wind  stream-function  values 
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Fig.  1 


1:  Contours  of  experimental  total  moisture  mixing  ratio 
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Fig.  12:  Contours  of  empirical  cloud  profiles 
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Fig.  13:  Sea  surface  sensible  and  latent  heat  flux  values 
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The  hybrid  model  developed  in  the  program  preserves  details  of  turbulence  of  the  second- 
order-closure  model  and  calculation  efficiency  of  the  large-eddy  simulation.  To  establish  confidence 
in  prediction  accuracy,  the  cloud  profiles  predicted  in  Section  4.3  using  the  hybrid  model  shown  in 
figure  14A  are  compared  with  the  empirical  dataset  presented  in  Section  5.1  plotted  in  figure  14B. 
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Fig.  14A:  Predicted  Cloud  Profiles 


Fig.  14B:  Empirical  Cloud  Profiles 


6.  CONCLUSION 


# 
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For  goals  of  research  at  UDC  to  identify  physical  processes  that  determine  the  dynamics  of 
marine  cloud  layers  and  to  quantify  roles  of  turbulence,  convection  and  thermal  radiation  that  play 
in  formation,  dissipation  and  stability  of  the  marine  cloud  layers,  the  author  and  his  cc-workers  at 
UDC  have  endeavored  to  achieve  the  objectives: 

•  To  advance  turbulence  models  using  efficient  numerical  schemes, 

•  To  make  available  computer  simulation  programs  for  predicting  stability 
of  marine  cloud  layers, 

•  To  retrieve  web-site  data  of  marine  cloud  layers  from  reliable  sources, 

•  To  establish  confidence  in  models  through  experimental  verification,  and 

•  To  yield  insights  into  the  cloud’s  physical  processes. 

For  these  objectives,  two  turbulent  models,  using  the  second-order-closure  and  the  LES 
model,  respectively,  have  been  formulated.  A  finite-volume  procedure  has  been  developed  to  solve 
the  resultant  equations.  In  addition,  to  retain  turbulence  details  of  the  second-order  turbulence  and 
calculation  efficiencies  of  the  LES  model,  a  hybrid  model  has  been  assembled.  It  uses  a  multi¬ 
dimensional  framework  of  the  LES  model  but  uses  a  set  of  unidirectional  second-order  equations  to 
calculate  the  turbulent  intensity  values. 

For  experimental  verification  of  the  hybrid  model,  theoretical  results  have  been  compared  with 
relevant  empirical  data  retrieved  from  reliable  sources.  Several  data  sources  such  as  TOGA-COARE 
(Tropical  Ocean  Global  Atmosphere  -  Coupled  Ocean  Atmosphere  Response  Experiment),  GCE 
(Goddard  Cumulus  Ensemble),  FIFE  (The  First  ISLCP  Field  Experiment),  GDAAC  (Goddard 
DAAC)  and  others  were  evaluated  based  on  the  time-interval,  grid-size  and  domain  criteria. 


29 


According  to  these  criteria,  two  sources  of  data  TOGA-COARE  (experimental)  and  GCE  model 
(semi-empirical)  were  examined  in  detail.  Through  cooperation  with  the  NASA  scientists  and  a 
permission  of  using  their  CRAY  computer,  several  sets  of  the  GCE  data  were  decoded  and 
transmitted  to  the  author.  These  data  have  been  used  to  verify  the  UDC’s  model  simulation  results. 
Figure  14A  shows  snapshots  of  cloud  profiles  that  have  been  predicted  by  the  hybrid  model  under  the 
same  wind-stream-function,  air/sea  interface  sensible  and  latent  heating  fluxes  and  other  boundary- 
condition  values  that  were  provided  by  the  empirical  data.  Although  predictions  are  not  expected  to 
have  a  complete  agreement  with  the  experimental  data  (as  precipitation  is  not  considered  in  the 
present  theory),  a  good  degree  of  similarity  can  be  observed  between  the  predicted  cloud  profiles 
shown  in  figure  14A  and  those  empirical  cloud  profiles  shown  in  figure  14B. 
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ABSTRACT 

A  study  has  been  made  to  identify  physical  processes  that 
determine  the  dynamics  of  the  marine  atmosphere  and  its  cloud  layers. 
The  effort  starts  from  formulation  of  the  governing  equations  for 
conservation  of  momentum,  enthalpy  and  moisture  of  atmospheric  air. 
The  turbulence  transport  equations  are  derived  to  calculate  the 
Reynolds-stress  and  turbulence-flux  correlations  that  appear  in  the 
conservation  equations.  The  thermal  radiation  equations  are  derived 
to  calculate  the  radiation  flux  divergence.  A  virtual  dry  potential 
temperature  value  is  introduced  to  account  for  the  effects  of  water 
vapor  condensation  and  cloud  formation  on  the  transport  of  the 
turbulence  quantities.  A  simulation  computer  program  is  developed  to 
calculate  dynamic  responses  of  the  marine  atmosphere  and  its  cloud 
layers.  Results  of  numerical  experiments  have  led  us  to  a  better 
understanding  of  the  effects  of  heat  and  moisture  exchange  between 
air  and  sea,  cloud-top  radiative  cooling,  and  interactions  of 
atmospheric  turbulence  and  thermal  radiation  on  the  dynamics  of 
marine  cloud  layers. 


CONSERVATION  AND  TURBULENCE  EQUATIONS 
Let  U,  V  and  W  be  the  mean  velocity  components  in  the  east,  north 
and  vertical  (x,y,z)  directions  and  at  the  time  t,  FR  the  radiation  flux, 
0  the  mean  moist-air  potential  temperature  defined  as 
[T4<gz+LQvyC|J,  and  Q  the  mean  total  moisture  mixing  ratio  of  vapor 
and  liquid  water  (Qv“KJi)-  The  lowercase  symbols  u,  v,  w,  0  and  (*> 
represent  the  turbulent  fluctuation  quantities  corresponding  to  their 
respective  mean  quantities  U,  V,  W,  0  and  Q.  In  addition,  the 
following  definitions  will  be  used:  L  for  the  water  latent  heat  of 
vaporization,  Cp  the  constant  pressure  specific  heat,  T  the  temperature, 
g  die  gravitational  acceleration,  f  the  Coriolis  parameter,  0V  the  virtual 
dry  potential  temperature  defined  as  [T(l+1 .609Qv-Q)+gz/C,J,  and  p 
the  density.  Far  above  a  sea  surface,  the  geostrophic  balance  is 
assumed.  It  gives  the  following  relationship  between  the  pressure  field 
and  horizontal  velocity  components  Ug  and  Vf  (Holton,  1972): 


(1) 


INTRODUCTION 

Advances  in  turbulence  modeling  and  computational  techniques 
have  made  it  a  common  practice  to  study  numerically  the  turbulent- 
flow  problems  (Meilor  and  Yamada,  1974;  Moeng  and  Arakawa, 
1980;  Moeng  and  Randall,  1984;  Chi,  1991).  Several  studies  of 
turbulent  flows,  e.g.,  using  the  eddy-viscosity  model  for  a  tomado-like 
vortex  (Chi  and  Jih,  1974),  the  large-eddy-simulation  model  for  end- 
wall  boundary  layers  of  intense  vortices  (Chi,  1977),  the  €-k  model  for 
vortex  flows  over  the  water  surface  (Chi,  1987)  and  the  second-order 
closure  model  for  heat  and  moisture  transfer  in  marine  atmospheres 
(Chi,  1994),  have  been  reported  earlier.  In  this  paper,  a  study  of  the 
simultaneous  effects  of  turbulent  heat  and  moisture  transfer,  water 
vapor  condensation  and  thermal  radiation  on  the  dynamics  of  marine 
cloud  layers  is  presented. 


i  dp 

P  dy 


(2) 


In  addition,  equilibrium  of  velocity  and  temperature  between  phases 
is  assumed,  and  a  Boussinesq  approximation  is  used.  Consequently, 
the  conservation  equations  for  momentum,  enthalpy  and  moisture  of 
atmospheric  air  can  be  written  as: 


dt  J  *  dz  dz 


(3) 
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dz  dz 


(4) 


dQ_  pdQ  dwffi,  1  3Fr 
dt  dz  dz  p Cp  dz 


Here,  E,  represents  the  emissi  vity  value  E  at  the  height  z  where  the  FR 
value  is  being  calculated. 

NUMERICAL  PROCEDURE  AND  COMPUTER  • 

ALGORITHM 

The  conservation  and  turbulence  transport  equations  discussed 
above  can  be  written  in  the  following  general  form: 


)+<b=0,  within  domain  (10)  - 

dt  jdXj  dXj  dxj  • 


dQ  __jj/dCl  _dwu> 
dt  dz  dz 


(6) 

and  their  boundary  and  initial  conditions  may  be  written  as: 


The  turbulent  transport  equations  (Chi,  1991)  derived  from  the 
second-order-closure  assumption  can  be  used  to  solve  a  set  of  the 
fifteen  turbulent  variables  P,  V1,  w7,  u"v ,  u  w,  v  w,  uF,  75,  wF, 
ug>,  vu>,  wo),  F7,  Fg>,  and  w7;  they  include  the  variables  wu,  wv, 
\IF^  and  wuin  Eq.  (1-4).  To  allow  for  water  vapor  condensation  and 
cloud  formation,  the  buoyancy  terms  associated  with  fluctuation  of  the 
virtual  dry  potential  temperature,  uF^  VFv,  wF  v  FF  „and  wF  v 
appear  in  flic  turbulence  transport  equations;  they  may  be  calculated 
from  the  predicted  turbulence  variables,  using  the  following  equations 
for  the  clear  and  cloudy  layers,  respectively: 

I5>se*(0.6094>-1)-^  (7) 


(8) 


Here  £  represents  any  turbulence  quantities,  u,  v,  w,  etc.;  <J>  is  equal  to 
CpTL;  y  is  equal  to  (UCJdOT/dT,  and  superscript  •  represents  values 
at  the  saturation  state. 


THERMAL  RADIATION  EQUATIONS 

The  radiation  calculation  requires  a  large  amount  of  time.  Simple 
emissi  vity  approach  using  the  gray-gas  assumption  has  been  used. 
The  theory  discussed  by  Sparrow  and  Cess  (1966)  was  followed.  By 
applying  the  boundary  conditions  of  the  blackbody  radiation  at  the  sea- 
surface  and  zero  downward  radiation  at  the  top  of  the  atmosphere,  and 
letting  E  be  the  mixed  water- vapor-and-cloud  emissi  vity  value,  the 
thermal  radiation  flux  FR  can  be  calculated  by  the  equation: 


FR=°T04e 


2 

2 


(9) 


y—+aa+b=0,  on  boundary  T 
dn 


(ID  • 


a(x1y7z,0)=a0(j:1yTz),  at  initial  time  zero  (12) 


where  a(z,t)  are  mean  values  of  velocity  components,  enthalpy,  0 
moisture  and  turbulent  correlations',  4>(z.t)  values  are  sources/ sinks  for 
a;  the  Four  items  in  Eq.  (10)  describe  the  system's  inertia,  convection, 
diffusion  and  source/sink  strength,  respectively;  and  values  of  a'  and 
*b'  in  Eq.  (1 1)  may  be  specified  to  enforce  the  appropriate  boundary 
conditions. 

A  finite-element  procedure  (Chi,  1994)  is  used  for  the  discretization 
of  spatial  variables.  The  flow  domain  to  be  simulated  is  divided  into  9 
small  finite  elements.  Multiplying  both  sides  of  Eq.  (10)  by  a 
weighing  function  G)  and  integrating  over  the  finite  element  Qc,  the 
weighted  residual  equation  set  can  be  derived: 


r  iVa 


&.  &. 


(13) 


Integrating  the  second-order  term  in  Eq.  (13)  by  parts  and  substituting 
Eq.  (1 1 )  for  normal  gradients  at  the  element's  boundary  result  in  the 
following  equation: 

f  [o*L.T^ii£.U0>i*. *)<*♦/  (a<ja*b<j>)da -0  (14) 

Jo  dt  dx  dx  '  dx  Tr 

•  /  /  t 

The  weighing  functions  g>  and  dependent  variables  cc(x,y,z,t)  are 
int-^-olated  in  each  element  Q*  as  follows: 


(16) 


m 
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where  i|ra  denotes  the  interpolation  functions  at  the  nth  node  of  each 
finite  element,  and  and  Q  „are  the  dement  nodal  values  of  the 
variable  a  and  the  weighing  function  <*>,  respectively. 

Substituting  Eq.  (15)  and  (16)  into  Eq.  (14)  yields  the  following 
finite-element  equation: 


C  — -+(K+U+H,JA+Ln=  0 

HOt  ^  '  fWt  IWf  MW  ft  ft 


where  C,^  K^,  and  are  the  mth-row  and  nth-column 
elements  of  square  matrices  which  describe  distributions  of  the  finite 
element's  capacitance,  diffusivity,  convectiveness  and  surface  flux, 
respectively,  and  L„  is  the  nth  element  of  a  load  vector  for  the  finite- 
element  source/sink  terms. 

Assembling  the  finite-element  Eq.  (17)  over  the  whole  flow  field 
y  Is  a  set  of  ordinary  differential  equations: 

[qff^+[K+r+fl]{A}+[I]=0  (18) 

at 

Where  the  vector  {A}  contains  all  node  point  variables  lying  within 
the  t  ow  domain  Q;  the  matrices  [C],  [K],  [U]  and  [H]  represent  global 
distributions  of  capacitance,  diffusiveness,  convection  and  surface 
conductance,  respectively;  and  [L]  is  the  global  load  vector. 

The  finite-element  procedure  described  above  has  trans-formed  the 
initial  value  partial  differential  conservation  and  turbulence  transport 
equations  into  a  large-order  system  of  the  ordinary  differential 
equation  set  (18).  Solution  of  (18)  may  start  from  a  Taylor  series 
expansion  of  {A}  about  time  t 

atA)  .  d(A) 

(A)  ,={A}+Af[6 - —-(1-6)-—=-]  (19) 

dt  at 

where  0  is  an  implicit  coefficient  Multiplying  both  sides  of  Eq.  (19) 
by  capacitance  matrix  [C],  Eq.  (19)  may  be  re-written  as: 

a  (A)  ,  dtAJ 

(F(A„.1)l=[a[(A)..l-(A,)-(A7)e— -(Arid  -9)—]  (20) 

With  expressions  fortune  derivatives  provided  by  Eq.  (18),  Eq.  (20) 
becomes: 

aU-iao 

+ Af  (1  -0)([Jf+  U +H]„{A))I +{Ll„)  (21) 

The  equation  set  (21)  may  be  solved  for  {A}^  values  at  time  V, 
from  the  known  { A}  a  values  at  time  ^  by  a  Newton’s  iterative  process. 
A  computer  program  has  been  written  to  facilitate  the  process.  The 
resultant  computer  program  can  be  used  to  simulate  the  dynamics  of 


the  marine  atmosphere  and  its  cloud  layers.  Presented  below  are  two 
examples  of  simulation  results. 

RESULTS  AND  DISCUSSION 

Firstly,  a  numerical  experiment  was  made  to  show  responses  of  a 
marine  atmosphere  to  variations  of  the  sea  surface  temperature.  Here 
the  geostrophic  wind  velocity  aligned  in  the  x  direction  is  set^at  a 
constant  value  of  10  m/s;  the  Coriolis  parameter  f  is  set  at  1  -Oxlp  s  * 
and  the  sea  surface  roughness  parameter  Zp  is  set  at  0.1  m.  Initially, 
file  vertical  profiles  of  the  potential  temperature  are  chosen  to  be  well- 
mixed  at 290  °K  and  have  a  total  mixing  ratio  value  corresponding  to 
100  percent  relative  humidity  at  the  top  of  the  marine  planetary 
boundary  layer  (MPBL).  Just  above  the  MPBL  top  at  z  equal  to  1  km, 
there  is  an  inversion  layer 400  m  thick,  in  which  temperature  increases 
and  mixing  ratio  decreases  at  0.02  °K/m  and  8x10  m  ,  respectively. 
Above  the  inversion  top,  the  gradient  of  temperature  is  -0.005  °K/m 
and  that  of  the  mixing  ratio  is  zero.  Subsequently,  the  sea-surfacc 
temperature  is  allowed  to  vary.  Numerical  simulations  are  made  to 
experiment  on  the  dynamic  responses  of  the  marine  atmosphere  to 
variations  of  the  sea-surface  temperature. 

Figure  1  shows  a  hypothetical  cyclic  sinusoidal  temperature 
variation  of  die  sea  surface.  Calculations  were  allowed  to  proceed  for 
five  days  with  the  sea-surface  temperature  repeating  cyclically  every 
twenty  four  hours.  Shown  in  Fig.  2  through  4  are  predicted  cyclic 
variations  over  a  twenty-four-hour  period  of  the  mean  wind-velocity 
and  air-potential-temperature  distributions  at  different  time.  Predicted 
responses  of  turbulent  energy,  horizontal  components  of  the  Reynolds 
stress  and  vertical  component  of  the  enthalpy  flux  are  shown  in  Fig. 

5  through  8. 

Many  interesting  characteristics  of  the  marine  planetary'  boundary 
layer  (MPBL)  can  be  observed  in  these  figures.  It  can  be  seen  in  Fig. 

1  that  the  sea-surface  temperature  was  set  to  rise  during  the  hours  of 
ax  to  eighteen  and  fall  during  the  hours  of  eighteen  to  six.  As  the  sea- 
surfacc  temperature  falls  (rises),  stability  of  the  boundary  layer 
increases  (decreases).  The  trends  can  be  seen  in  Fig.  2  and  3  by 
decreasing  in  the  marine  surface  boundary  layer  thickness  during  the 
hours  of  eighteen  to  six  and  increasing  in  the  marine  surface  boundary 
layer  thickness  during  the  hours  of  eighteen  to  six.  In  particular, 
instability  of  die  marine  thermal  boundary  layer  can  be  observed  in  the 
vicinity  of  the  eighteenth  hour  when  the  sea-surface  temperature  is  at 
its  maximum.  It  can  be  observed  also  in  Fig.  5  where  contours  of  the 
turbulent  energy  values  at  different  height  and  time  are  shown. 
Instability  of  die  MPBL  in  die  vicinity  of  the  eighteenth  hour  can  be 

seen  also  in  this  figure.  Although  it  can  be  observed  in  this  figure,  the 
maximum  thickness  of  die  MPBL  in  the  experiment  does  not  occur  at 
the  exact  time  of  the  maximum  sea-surface  temperature  at  the 
eighteenth  hour  (but  at  approximately  the  21st  hour)  because  of  delay 
in  the  responses.  Similarly,  the  effects  of  stability  on  distribution  of 
the  horizontal  components  of  the  Reynolds  stress  and  the  vertical 
component  of  the  heat  flux  can  be  seen  in  Fig.  6  to  8.  Again, 
instabilities  of  the  MPBL  can  be  seen  in  these  figures  by  observing  the 
maximum  stress  and  flux  values  and  increasing  in  the  MPBL 
thickness  in  the  vicinity  of  the  eighteenth  hour. 

For  the  second  case  of  numerical  experiments,  its  initial  conditions 
at  the  zeroth  hour  are  obtained  from  the  eighteenth-hour  results  of  the 
first-case  expieriments  presented  above.  In  this  case,  the  sea-surfacc 
temperature  is  changed  at  die  zeroth  hour  to  298*K,  and  it  is 
maintained,  thereafter,  at  this  constant  value.  The  sea-surface-air 
moisture  mixing  ratio  is  maintained  at  the  corresponding  saturated 
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Fig.  3  Y-Component  mean  wind  velocity,  V  m/s  Fig.  4  Mean  air  potential  temperature,  0  *K 


TTm«,  hours 


Fig.  9  Contours  of  predicted  cloud  intensity  of 
the  liq  -  id-water  mixing  ratio  in  g/kg 


state.  Plotted  in  Fig.  (9)  are  contours  of  the  predicted  cloud  intensity 
of  the  liquid-water  mixing  ratio  values  in  g/kg.  It  can  be  seen  in  this 
figure  that  the  cloud  starts  to  form  at  the  14th  hour.  Its  depth  grows 
firstly  downward  because  of  the  upward  transport  of  the  water  vapor 
from  the  sea  surface  and  the  upward  drop  in  the  static  temperature  of 
the  air.  From  the  25th  hour  onwards,  because  of  the  cloud-top 
radiative  cooling,  the  depth  of  the  cloud  increases  at  the  top.  Also,  it 
can  be  seen  in  this  figure  that  the  steady  state  of  the  cloud  layer  was 
reached  at  about  the  48th  hour. 


CONCLUSION 

In  conclusion  a  numerical  study  of  turbulence  in  marine 
atmospheres  has  been  made.  Results  of  predictions  have  led  us  to  a 
better  understanding  of  the  effects  of  air/sca  heat  and  moisture 
exchange,  cloud-top  radiative  cooling,  and  interactions  of  atmospheric 
turbulence  and  thermal  radiation  on  the  dynamics  ot  manne  cloud 
layers. 
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Goals  of  our  atmospheric  research  at  UDC  are  to  identify  physical 
processes  that  determine  the  dynamics  of  the  cloud  layers  and  to 
quantify  the  roles  of  turbulence,  convection  and  thermal  radiation 
that  play  in  formation,  dissipation  and  stability  of  the  cloud  layers. 
Our  immediate  objectives  are  to  advance  theoretical  models,  use 
efficient  numerical  schemes  and  develop  computer  programs  to 
simulate  the  marine  cloud  layers.  Comparison  of  computer  results 
with  published  field  observations  will  yield  insights  mto  the  cloud- 
layers’  physical  processes.  While  a  companion  paper  FEDSM98- 
4954  develops  a  large-eddy  turbulent  model  to  simulate  stability  of 
the  marine  cloud  layers,  this  paper  describes  a  second-order  closure 
turbulent  model  to  simulate  turbulent  mixing  of  the  moist 
atmospheric  air  in  general  and  the  marine  atmospheres  in  particular 

introduction 

Turbulent  mixing  of  heat,  moisture  and  momentum  plays  a 
dominant  role  in  atmospheric  processes  that  are  of  interests  to 
scientists  of  different  disciplines:  physicists,  meteorologists,  and 
environmentalists.  Advances  in  turbulent  modelmg  will  lead  to 
better  understanding  of  climatic  phenonmena1A*  -  cloud  stability, 
tornado  motion,  and  severe-storm  formation  -  and  environmental 
qualities’’10,13  -  air  pollution,  acid  deposition  and  global  warming. 
Consequently,  improvements  can  be  made  in  climatic  prediction  and 
environmental  control. 

Much  of  our  understanding  of  mixing  processes  in  atmospheres 
has  come  from  careful  observation  and  sound  theoretical  modelmg. 
While  it  is  yet  impossible  to  have  a  generalized  turbulence  theory 
for  universal  phenomena,  semi-empirical  models  with  different 
degrees  of  complexity  have  been  developed  with  confidence  in 
simulating  numerous  practical  phenomena:  the  mixing-length  eory 


has  been  used  to  model  turbulent  boundary  layers  on  flat  plates," 
the  second-order  diffusion  model  has  been  used  for  simulating  the 
planetary  boundary  layers,6  and  large-eddy  models  have  been  used 
to  study  the  meso-scale  turbulence  in  atmospheres.  This  author 
has  used  the  eddy-viscosity  model  to  simulate  a  tomado-like  vortex^ 
the  large-eddy-turbulence  model  for  end-wall  boundary  layers  of 
intense  vortices,  the  €-k  model  for  vortex  flow  over  the  water 
surface,  and  the  second-order  closure  model  for  the  marine  cloud 
layer.  It  is  intended  in  this  study  to  present  a  hybrid  treatment  for  the 
atmospheric  turbulence  that  uses  a  second-order  closure  turbulent 
diffusion  model  for  simulating  atmospheric  mixing  layers  and  a 
large-eddy  turbulent  model  for  simulating  meso-scale  entrainment 
in  the  marine  atmosphere.  Presented  in  this  paper  is  a  second-order 
diffusion  model  for  simulating  moisture  mixing  in  the  marine 
atmosphere;  a  large-eddy  model  for  simulating  the  cloud-layer 
entrainment  will  be  presented  in  a  companion  paper. 

NOMENCLATURE 

C  =  specific  heat 

A^B  C  D  E  &  F  values  in  turbulent  flux  equation 

’=  second-order  closure  equations’ empirical  constants 

f  =  geostrophic  Coriolis  coefficient 

A.  R  =  thermal  radiation  flux 

g  =  gravitational  acceleration 

L  =  water-vapor  latent  heat  of  vaporization 

q  =  square-root  value  of  squared  mean  of  the  turbulent 

fluctuating  velocity 
t  =  time 

T  =  static  temperature  of  atmospheric  air 

(u,v,w)  =  turbulent  fluctuating  velocities  in  (x,y,z)  directions 
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(U,  V,W)  =mean  wind  velocity  components  in  (x,y,z)  directions 
(Ug,Vg)  =geostrophic  wind  components  in  (x,y)  directions 
(x,y,z)  =east-,  north-  and  vertical-direction  coordinates 
a  =a  generalized  variable  in  a  numerical  scheme 
P  -  buoyant  coefficient 

0  =  turbulent  fluctuating  value  of  8 

0V  =  turbulent  fluctuating  value  of  8V 
8  =  mean  potential  temperature  =  T-Kgz+LQJ/Cj, 

©v  =  mean  virtual  potential  temperature 
=  T( 1+ 1 . 60 9£lv-Q)+gz/Cp 
K  =  Prandtl  mixing  length  constant 

A  =  turbulent  length  scale 

(i>  =  turbulent  fluctuating  value  of  Q 

Cl  =  mean  total  moisture  ratio 

Cl,  =  mean  liquid  moisture  ratio 

Q*  =  mean  vapor  moisture  ratio 

Overlined  values 

=  turbulent-flux  values 
Other  symbols  in  numerical  procedures 
=  those  defined  in  the  text 

CONSERVATION  AND  TURBULENCE  EQUATIONS 

When  these  assumptions  are  made:  (1)  at  far  above  a  sea  surface, 
the  geostrophic  balance  is  maintained,  (2)  velocity,  temperature  of 
vapor  and  liquid  moisture  are  in  equilibrium,  and  (3)  Boussinesq 
approximation  is  used,  the  conservation  equations  for  momentum, 
enthalpy  and  total  moisture  of  atmospheric  air  can  be  written  as: 


dU -fry ~y  \  (1) 

dt  1  dz  dz 

W._AUU)_d™_wdV  (2) 

dt  g  dz  dz 


de_wdedwd  1  dFR 

dt  dz  dz  pC  dz 

p 


dCl  jpdCl  0Wi) 
dt  dz  dz 


(4) 


_  _  • 

(8) 

dt  dz  dz  dz  A 

(9) 

dt  dz  dz  dz  A  W 


dt  dz  dz  dz  dz  X 

(10) 

0v0  d  /y4  .  0v0.  — 08  -—ndV  ^q-s 

- -—{AqX - )~vw — -w0— 

dt  dz  dz  dz  dz  A 

(11) 

dt  dz  dz  dz  A 

(12) 

^L=l{AqX— )-2*0— -FI# 

dt  dz  dz  dz  X 

(13) 

dut 0  d/A  .01/(1).  — 0Q  — dU  nq  — 

- = — (AqX - )-uw — “>vu>— — E-f-i/o) 

0/  dz  dz  dz  dz  A 

(14) 

0vg)  0..  ,  0vo).  — 0Q  —dV  vq— 

- ---(AqX - )~vw — -Hti)-— -jK-7-vo) 

dt  dz  dz  dz  dz  A 

(15) 

^--±(2AqX^)-^B^g^-Elw 
dt  dz  dz  dz  A 

(16) 

^--^AqX^-ylwQ-Flu1 

dt  dz  dz  dz  A 

(17) 

00( 1)  6/1  «  00(0 .  00  a  0fl  rr  q  n 

- =^-(AqX - )-nti) - w0 — -F—  00) 

dt  dz  dz  '  dz  dz  X 

(18) 

In  above  equations,  values  for  overlined  turbulent  -flux  values  can 
be  calculated  by  the  turbulent  transport  equations  using  the 
assumptions  of  Mellor  and  Yamada’s  second-order  closure  model:7 


dtS 

dt 


dzx  H  dz  dz  r  3  '  3A 


(5) 


dt  dz  dz  dz  X  3  3A 


(7) 


In  above  equations,  A  is  a  characteristic  length  scale  that  is  equal  to 
the  value  of  the  Blackadar's  or  the  diffusion-length  scale  -  AB  or  AD  - 
whichever  is  the  smallest. 

In  equations  5  through  18,  the  time-derivative  terms  on  the  left- 
hand  Ji  -lj  model  the  transient  variation  of  turbulence  correlations.  £ 
The  second-order  derivative  terms  on  the  right-hand  side  model 
turbulent  diffusion.  While  the  production  of  turbulence  due  to 
buoyancy  is  modeled  by  terms  with  buoyant  coefficient  p,  the 
production  of  turbulence  due  to  friction  is  modeled  by  products  of 
second-order  correlations  and  gradients  of  mean  variables.  Terms 
with  coefficients  B,  C,  E  and  F  represent  the  turbulent  redistribution,  f 
Turbulent  dissipation  is  modeled  by  terms  with  coefficient  D. 
Coefficients  k.  A,  B,  C,  D,  E  and  F  have  been  determined  semi- 
empirically,3  having  the  values  equal  to  0.35,  0.21,  0.46,  0.053, 
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0.132,  0.44,  and  0.23,  respectively. 

NUMERICAL  PROCEDURE  AND  COMPUTER 
ALGORITHM 

The  conservation  and  turbulence  transport  equations  1  to  18 
described  above  can  be  written  in  the  following  general  form: 

§l+W^--—(y—)+6=0,  Within  domain  (19) 

dt  dz  dz  1  dz 

The  boundary  and  initial  conditions  may  be  written  as: 

y—+aa+b=09  On  boundary  V  (20) 

dn 

a(z,0)=a0(z),  At  initial  time  zero  (21) 


where  0  is  an  implicit  coefficient.  Multiplying  both  sides  of 
equation  24  by  capacitance  [C]  and  using  expressions  for  time 
derivatives  provided  by  equation  23,  equation  24  becomes: 

IFHQ U^-rlAU 

+A/0([^+U+/f]„+1<A}B+1  +<&U)  (25) 

+A/(1  -Q)([K+U+H]rtWn+{L)n) 

The  equation  set  25  may  be  solved  for  the  vector  {A}^  at  time 
t^,  from  the  known  (A}n  values  at  time  t  by  a  Newton's  iterative 
process,  and  a  computer  program  has  been  written  to  facilitate  the 
process.  The  resultant  computer  program  can  be  used  to  simulate 
mixing  processes  in  the  marine  atmosphere.  Presented  below  is  an 
example  of  computer-simulation  results. 


In  above  equations,  a(z,t)  values  are  mean  values  of  velocity 
components,  enthalpy,  moisture  and  turbulent  correlations;  4>(z,t) 
values  are  sources/sinks  for  a.  The  Four  items  in  equation  19 
describe  the  system’s  inertia,  convection,  diffusion  and  source/sink 
strength,  respectively,  'a'  and  V  values  in  equation  20  are  used  to 
define  the  appropriate  boundary  conditions. 

A  finite-element  procedure4  was  used  for  descritization  of  spatial 
variables.  The  flow  domain  to  be  simulated  was  divided  into  small 
finite  elements.  Multiplying  both  sides  of  equation  1 9  by  a  weighing 
factor  and  integrating  over  the  finite  element,  a  weighted  residue 
equation  set  can  be  derived.  Then,  substituting  into  the  residue 
equation  a  set  of  appropriate  interpolation  functions  for  the  weighing 
factor  and  the  element’s  variables  resulted  the  following  finite- 
element  equation  set: 

C  UK  +U  +H  )A  +L  =0  (22) 

“"l/T  mn  mn  mn  n  n 


Where  Cm,  K^,  li,  and  ft  are  the  mth-row  and  nth-column 
elements  of  square  matrices  which  describe  distributions  of  the  finite 
element's  capacitance,  difiusivity,  convectiveness  and  surface  flux, 
respectively,  and  is  the  nth  element  of  a  load  vector  for  the  finite- 
element  source/sink  terms.  Finally,  assembling  the  finite -element 
equation  22  over  the  whole  flow  field  yielded  a  set  of  ordinary 
differential  equations: 


U+H]{ A) + [L]  =0 
dt 


(23) 


Where  the  vector  {A}  contains  node-points  variables  lying  within 
the  flow  domain  Q,  [L]  is  the  global  load  vector,  and  the  matrices 
[C],  [K],  [U]  and  [H]  represent  global  distributions  of  capacitance, 
diffusiveness,  convection  and  surface  conductance,  respectively. 

The  finite-element  procedure  described  above  has  transformed  the 
initial  boundary  value  partial  differential  equations  1  through  1 8  into 
a  large-order  system  of  the  ordinary  differential  equation  set  23. 
Solution  of  23  may  start  from  a  Taylor  series  expansion  of  {A} 
about  time  t^: 

3(A)  .  3(A) 

<AU  =lA>M0-^+(l-0)-gf]  (24) 


RESULTS  AND  DISCUSSION 

The  above-described  computer  program  has  been  used  to  predict 
transient  exchange  of  moisture  in  the  sea/air  interface,  mixing  of 
moisture  in  the  marine  atmosphere,  and  formation  of  the  cloud  layer. 
Graphs  plotted  in  figure  1  were  reported  in  a  previous  paper,  they 
were  used  as  initial  conditions  for  this  study. 


Fig.  1 :  Initial  Velocity-Components  and  Potential-Temperature 
Values 
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The  initial  conditions  shown  in  figure  1  were  obtained  on  the 
assumptions  that  the  geostrophic  wind  in  the  x-direction  was  at 
lQm/s,  and  Coriolis  parameter  f  is  equal  to  1 .0x1  O*4  s*1,  the  potential 
temperature  was  chosen  to  be  well-mixed  at  290K,  and  the  total 
moisture  mixing  ratio  corresponding  to  100  percent  relative 
humidity  at  top  of  the  marine  planetary  boundary  layer  (MPBL). 
Just  above  the  MPBL  top  at  height  equal  to  1  km,  there  was  an 
inversion  layer  of 400-m  thick,  in  which  temperature  was  increasing 
at  0.2  K/m  and  moisture  mixing  ratio  was  decreasing  at  0.008 
g/kg.m,  respectively.  The  sea-surface  temperature  was  allowed  to 
vary  diumally  within  the  range  of 280  to  290F.  Conditions  shown 
in  figure  1  were  for  the  instance  when  the  sea  surface  was  at  290  KL 
For  the  simulation  run  using  initial  conditions  shown  in  figure  1 , 
it  was  assumed  that  the  sea  surface  temperature  was  raised  abruptly 
to  and  then  maintained  at  298  K,  and  the  air  total  moisture  ratio  at 
the  sea  surface  was  maintained  at  the  saturation  state.  After  a 
simulation  run  over  a  period  of  forty  hours,  a  large  body  of  data  was 
generated.  Plotted  in  figures  2  are  predicted  contours  of  temporal 
mean  physical-property  values  of  air  in  the  marine  atmosphere; 
plotted  in  figures  3  are  predicted  contours  of  Reynolds-stress  and 
turbulent-flux  values  in  the  marine  atmosphere.  Many  interesting 
characteristics  of  mixing  in  the  simulation  domain  can  be  observed 
in  these  contours.  Firstly,  rapid  interaction  at  the  air/sea  interface 
can  be  observed  during  the  initial  period  of  zero  to  600  minutes.  It 
is  followed  by  a  calmer  development  of  the  marine  planetary 
boundary  layer  for  about  600  minutes.  Then,  during  the  next  600 
minutes  transfer  of  enthalpy  and  moisture  continues,  as  can  be 
observed  from  the  predicted  contours  of  the  turbulent  thermal-  and 
moisture-flux  values  shown  in  figures  3C  and  3D.  Also,  can  be 
observed  in  Fig.  2D  is  the  formation  of  cloud  starting  at  the  200th 
minute,  rapid  deepening  of  the  cloud  layer  during  the  second  interval 
of 200  to  1200  minutes,  and  slower  growth  of  the  cloud  layer  during 
the  period  of  1200  to  1800  minutes.  Finally,  a  steady  state  is 
established  at  around  the  2400th  minute. 

Plotted  in  figure  4  are  the  predicted  steady-state  conditions  at  the 
end  of  the  simulation  period  of  the  40th  hour.  These  graphs  shown 
in  figure  4  can  be  compared  with  those  in  figure  1  for  the  initial 
conditions.  Changes  that  have  been  made  during  these  forty  hours 
can  be  observed.  Firstly,  warm  and  moist  air  at  the  sea  surface  has 
resulted  in  an  unstable  mixing  layer.  It  can  be  observed  by 
comparison  of  the  mean  horizontal  velocity  components  plotted  in 
figures  1 A  and  4  A  of  the  thickening  of  the  unstable  boundary  layer 
in  figure  4 A  due  to  the  intensed  turbulent  exchange.  It  can  also  be 
observed  in  figure  4H  of  the  increased  turbulent  kinetic  energy  in 
comparison  with  the  corresponding  turbulent-energy  values  plotted 
in  figure  1.  In  addition,  variations  of  profiles  for  the  mean  wind 
velocity,  mean  potential  temperature,  mean  total-moisture  ratio, 
mean  liquid-moisture  ratios,  Reynolds-stress,  turbulent  thermal  flux, 
turbulent  moisture  flux  and  turbulent  kinetic  energy  can  be  observed 
in  figure  4. 

CONCLUSION 

In  conclusion,  a  numerical  study  of  turbulent  mixing  of  moist  air 
in  the  marine  atmosphere  has  been  made.  Results  of  predictions 
have  led  us  to  a  better  understanding  of  the  turbulent  mixing  process 
in  the  atmosphere.  However,  complexity  of  the  turbulent  model  has 
limited  our  study  to  a  one-dimensional  model.  Physics  learned  in 
this  study  has  led  us  to  the  development  of  a  multidimensional  large- 


Time,  minute* 


Fig.  2:  Predicted  Contours  of  Mean  Atmospheric  Quantities 
2A  -  Mean  Wind  Velocity,  m/s 
2B  -  Mean  Potential  Temperature,  K 
2C  -  Mean  total  moisture  mixing  ratio,  g/kg 
2D  -  Mean  liquid  moisture  mixing  ratio,  g/kg 
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Fig.  3:  Predicted  Contours  of  Turbulent  Fluxes 
3A  -  Turbulent  kinetic  energy,  mA2/sA2 
3B  -  Principal  Reynolds  stress,  mA2/sA2 
3C  -  Turbulent  thermal  flux,  m.K/s 
3D  -  Turbulent  moisture  flux,  m.kg/kg.s3 


Fig.  4:  Predicted  Mean  Quantities  and  Turbulent  Fluxes 
4A  -  Mean  velocity;  4B  -  Reynolds  stress; 

4C  -  Potential  temperature;  4D  -  Thermal  flux; 

4E  -  Total  moisture  ratio;  4F  -  Moisture  flux; 

4G  -  Liquid  moisture  ratio;  4H  -  Turbulent  energy. 
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eddy  turbulent  model  for  investigating  stability  of  the  marine  cloud 
layer.3 
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ABSTRACT 

Goals  of  our  atmospheric  research  at  University  of  the  District 
of  Columbia  (UDC)  are  to  identify  physical  processes  that 
determine  the  dynamics  of  the  cloud  layers  and  to  quantify  the  roles 
of  turbulence,  convection  and  thermal  radiation  that  play  in 
formation,  dissipation  and  stability  of  the  cloud  layers.  Our 
immediate  objectives  are  to  advance  theoretical  models,  use 
efficient  numerical  schemes  and  develop  computer  programs  to 
simulate  the  marine  cloud  layers.  Comparison  of  computer  results 
with  published  field  observations  will  yield  insights  into  the  cloud- 
layers’  physical  processes.  While  a  companion  paper  FEDSM98- 
4809  develops  a  second-order  closure  model  for  simulating 
turbulent  mixing  of  the  moist  marine  atmosphere,  this  paper  deals 
with  a  large-eddy  turbulent  model  for  simulating  stability  of  the 
marine  cloud  layers. 

INTRODUCTION 

Much  of  our  understanding  of  turbulent  processes  has  come  from 
careful  observation  and  sound  theoretical  modeling.  While  it  is  yet 
impossible  to  have  a  generalized  turbulence  theory  for  universal 
phenomena,  semi-empirical  models  with  different  degrees  of 
complexity  have  been  developed  with  confidence  in  simulating 
numerous  practical  phenomena:  the  mixing-length  theory  has  been 
used  to  model  turbulent  boundary  layers  on  flat  plates.11  the 
second-order  diffusion  model  has  been  used  for  simulating  the 
planetary  boundary  layers,7’8-9  and  large-eddy  models  have  been 
used  to  studythe  meso- scale  turbulence  in  atmospheres.10,12  This 
author  has  used  the  eddy-viscosity  model  to  simulate  a  tornado-like 
vortex,1  the  large-eddy-turbulence  model  for  end-wall  boundary 


layers  of  intense  vortices,2  the  e-k  model  for  vortex  flow  over  die 
water  surface,3,4  and  the  second-order  closure  model  for  the  marine 
cloud  layers.5*  This  study  of  the  author  uses  a  hybrid  treatment  for 
the  atmospheric  turbulence  that  employs  a  second-order  closure 
turbulent  diffusion  model  for  simulating  atmospheric  mixing  layers 
and  a  large-eddy  turbulent  model  for  simulating  convective 
entrainment  at  the  cloud  top.  The  stratus  cloud  plays  an  important 
role  in  climatic  dynamics;  it  has  stimulated  extensive  research. 
Second-  and  higher-order  turbulent  models  have  succeeded  in 
advancing  theoretical  understanding  of  the  marine  cloud  layers. 
Complexity  of  those  models  often  makes  long-term  simulation  of 
the  cloud  layer  over  an  extensive  period  prohibitively  expensive. 
A  hybrid  model  using  second-order  bottom-up  mixing  and  large- 
eddy  top-down  convection  will  improve  in  computational 
efficiency. 

NOMENCLATURE 

Symbols: 

a,b  =  coefficients  in  boundary-layer  or  numerical  equations 

C  =  specific  heat 

D  =  diffusion  parameter 

E  -  turbulent  kinetic  energy 

F  =  convection  flux  parameter 

g  =  gravitational  acceleration 

K  -  eddy  coefficient  values 

L  —  water- vapor  latent  heat  of  vaporization 

p  =  pressure 

t  =  time 

T  =  resolved  mean  static  temperature  of  atmospheric  air 
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(U,  V,W)  =  resolved  mean  wind  velocities  in  (x,y,z)  directions 
(u,v,w)  =  turbulent  fluctuating  velocities  in  (x,y,z)  directions 
(x,y,z)  -  east-,  north-  and  vertical-direction  coordinates 

a  =  a  generalized  variable  in  a  numerical  scheme 
P  =  buoyant  coefficient 

r  =  diflusivity  value  in  the  generalized  differential  equation 

y  -  difiiisivity  value  in  the  generalized  boundary  condition 

6  =  a  length  scale  defined  by  equation  1 8 

(Ax,Ay,A z)  =  node  points’  distance  in  (x,y,z)  directions 
€  =  a  turbulent  energy  dissipation  rate 

X]  =  a  moisture  parameter  defined  as  QJC?)dQ*/dT 
0  =  resolved  mean  potential  temperature  =  T-Kgz+LQJ/Cp 

0V  =  resolved  mean  virtual  potential  temperature 

=  T(l+1.609CVG)+gz/Cp 

Q  =  resolved  mean  total  moisture  mixing  ratio 

=  Q.  +  0, 

Q*  =  saturation  moisture  mixing  ratio  value  at  air  temperature 

A  =  a  defined  domain  under  simulation 

A.  =  turbulence-length  scales 

£  =  a  moisture  parameter  defined  as  CpT/L 

II  =  a  defined  domain  boundary 

p  —  air  density 

o  =  turbulent  Prandtl  number 

<f)  =  source  terms  in  generalized  differential  equations 

Suffices: 

app  =  for  intermediate  approximate  value  in  an  iterative  process 
e  =  for  turbulent  kinetic  energy 

t  =  for  liquid  water 

m  =  for  momentum 

o  =  for  a  reference  state 

(P^,WJB,D  “ 

for  a  variable  node  point  in  a  domain  element’s  center 
and  its  neighboring  node  pointes  at  the  east,  west,  top 
and  bottom 
(p,e,w,b,t)  = 

for  a  velocity  node  point  on  an  edge  of  the  domain  element 
and  its  neighboring  velocity  node  pointes  at  the  east,  west, 
top  and  bottom 
v  =  for  water  vapor 

0  =  for  enthalpy 

a)  =  for  moisture 


—=—[/:  (—+—)] +2-^/:m—)+pg(e-eo) 
dt  dxl  mKdz  a*  dz  mdz’ Vh  0 

vdw  wdw  1  dP 
dx  dz  po  dz 


(2) 


)+A(r  —  )-U— -F— 
dt  dx  0  ax  dz  0  dz  dx  dz 


(3) 


HQ.=Ajk  —)-U—-W—  (4) 

dt  a?  • dx  ’  dz  udz  ’  dx  dz 


dt  dx  dz  dx  E  dx  dz  dz  0O  dz 

+2K  {—?+K  (—+—f+2KJ—?]~e 
m  dx  "  dz  dx  m  dz 


(7) 


dJL+W-  0 

dx  dz 


(6) 


In  above  equations,  the  buoyancy  terms  associated  with  the  virtual 
dry  potential  temperature  have  been  defined  as  follows: 


0v=7’(l  +  1.6O9Qv-fi)+^-  (7) 

P 


dz  dz  C  dz 

p 

For  Clear  Air  Layers 


(8)  • 


CONSERVATION  AND  TURBULENCE  EQUATIONS 

When  these  assumptions  are  made:  (1)  the  Coriolis  force  is 
negligible,  (2)  velocity  and  temperature  of  vapor  and  liquid  moisture 
are  in  equilibrium,  and  (3)  Boussinesq  approximations  are  used,  the 
conservation  equations  for  momentum,  enthalpy  and  total  moisture 
of  atmospheric  air  can  be  written  as: 


dU  ,  d  ,K  dU.  3r„  ,dU  dW., 
dt  dx  dx  dz  dz  dx 

vdU  wdU  1  dP 
dx  dz  pedz 


(1) 


de_it-i. 609friae  ,L_dQ 
dz  1+q  dz  Cf  dz 
For  Cloudy  Air  Layers 


(9)  • 


The  rate  of  dissipation  within  the  grid  volume  and  the  subgrid  eddy 
coefficient  may  be  parameterized  through 


€=0.19 Em!X 


(10) 


Arm=0.58£1/2Jl 


(ID 
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K  2X 

^a=— (1+4— ) 

°e  *-s 


K  2X 

K  = — — (1  +— ) 
a  Ae 


source/sir.  ^  6ins  for  a.  Items  in  equation  19  describe  the 

(12)  variables’  time  derivative,  convection,  diffusion  and  source/sink 
strength,  respectively,  'a'  and  V  values  in  equation  20  are  used  to 
define  the  appropriate  boundary  conditions. 

A  finite-volume  difference  scheme13  was  used  for  descritization  of 

(13)  variables.  The  flow  domain  to  be  simulated  was  divided  into  small 
rectangular  elements;  shown  in  figure  1  are  examples  of  several 
such  elements.  It  can  be  seen  in  figure  1  that  variable  values  at  node 
points  of  the  elements  are  strategically  located.  The  node  point  for 

(14)  velocity  value  is  located  in  the  middle  of  the  rectangular  edge  and 
the  node  point  for  other  variables  is  located  in  the  center  of  the 
rectangular  box. 


where  A  is  a  minimum  of  the  Blackadar’s  length  A0,  difiusion  length 
Ad,  resolvable  length  scale  Xs: 


x  _  0.3525 

B  8+3.5 z 

(15) 

V-0.75O'» 

(16) 

Xs=(AxAxAz),rs 

(17) 

f~Emzdz 


5=  0 

(18) 

(-Emck 

And  the  terms  o0,  o^ ,  and  ,  which  are  the  turbulent  Prandtl 
numbers  of  enthalpy,  moisture  and  turbulent-energy  transports,  are 
equal  to  0.75, 0.75  and  1,  respectively. 


NUMERICAL  PROCEDURE  AND  COMPUTER 
ALGORITHM 

The  conservation  and  turbulence  transport  equations  1  to  5 
described  above  can  be  written  in  the  following  general  form: 


§1 + + fvil =. -  JL(A -i-(r^) +$> 

dt  dx  dz  dx  dx  dz  dz 
within  Domain  A 


(19) 


The  boundary  and  initial  conditions  may  be  written  as: 

y—+aa+b=0,  On  boundary  II  (20) 

dn 

a(xjj=0)=a0(xj)y  At  initial  time  zero  (21) 

In  above  equations,  a(x*z,t)  values  are  variable  values  of  velocity 
components.  U  and  W,  potential  temperature  0,  moisture  mixing 
ratio  Q,  and  turbulent  kinetic  energy  ‘E’;  4>(x>z,t)  values  are  the 


Fig.  1 :  An  Element  with  Neighboring  Variables  Shown 

Using  the  velocity  values  and  values  of  a  shown  in  the  figure  and 
employing  an  upwind  logic,  the  conservation  law  may  be  applied  to 
obtain  an  expression  for  the  value  of  a  at  the  node  point  P  in  terms 
of  the  a-values  at  its  neighboring  node  points: 


apap=aEaE+affff.w+a7aT+aBaB+b  (22) 


where  coefficient  values  ‘a*  and  ‘b’  may  be  expressed  in  terms  of 
the  difiusion  and  flux  parameters  defined  as  follows: 


Kte 

aE +Max[  -(£/ Az),0] 
i*xt 


(23) 


aw=^-+Mcoc[(UJ>z),Q] 


(25) 


KAx 

aT=—j~~  +Max{-(IJ 


(26) 


K.Ax 

aB=-±—+Max[(UbAx)fi] 


(27) 


AxAz 

aP=aE+aW+aT+aB+-£- 


(28) 
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A=<()fAxAz 


(29) 


half,  one,  five  and  ten  hours  from  the  start  of  the  simulation  run), 
dissipation  of  the  cloud  layer  can  be  observed  0 


In  addition,  in  solving  U  and  W  values  using  equation  22,  the 
initial  pressure  values  will  have  to  be  the  estimated  values; 
consequently,  the  equation  22  will  yield  initially  approximate 
and  Ww  values.  Improved  U  and  W  values  may  be  calculated  from 
their  approximate  values  using  equations: 


U  =U  +—  (/>-/>) 


U  =U  +  (P  ~PD) 

w  w,app  a  o  V  w  p 
pr0 


CONCLUSION 

In  a  companion  paper,6  a  second-order  closure  model  was 
developed  to  simulate  sea-to-air  heat  and  moisture  transfer  and  to 
observe  theoretically  the  growth  of  the  marine  cloud  layers.  Results 
of  predictions  have  led  us  to  a  better  understanding  of  the  turbulent  ® 
mixing  process  in  the  atmosphere.  However,  complexity  of  the  high- 
order  closure  model  had  limited  our  study  to  the  mixing  layer 
problems  in  one  dimension.  Physics  learned  in  that  study  has  led  us 
to  the  development  of  a  multidimensional  large-eddy  turbulent 
model  presented  in  this  paper.  The  extended  theory  has  enabled  us 
to  investigate  multidimensional  effects  of  convective  entrainment  at  # 
the  cloud  top  on  the  dynamics  of  the  marine  cloud  layers.  These 
quantitative  results  have  provided  us  at  UDC  with  theoretical  tools 
to  identify  processes  observed  physically  in  the  marine  cloud  layers. 


wrwim*~£~^pr~PT>  <32> 


Substituting  U  and  W  values  in  equations  30-33  into  the  continuity 
equation  6  yields  a  set  of  linear  equations  for  pressure  values  at  the 
solution  domain’s  node  points;  they  may  be  used  to  solve  for 
improved  pressure  values.  So  the  process  may  be  repeated  to  iterate 
alternatively  for  improved  values  of  velocity,  pressure  and  other 
variables  at  the  node  points. 

Using  the  linear  equations  discussed  above  for  variable  values  U, 
W,  9,  Q,  E  and  P  at  node  points,  a  computer  program  has  been 
written  to  simulate  dynamics  of  the  marine  cloud  layers. 

RESULTS  AND  DISCUSSION 

The  above-described  computer  program  may  be  used  to  study 
stability  of  the  marine  cloud  layers.  Graphs  plotted  in  figure  2  show 
characteristics  of  a  horizontally  uniform  cloud  layer  predicted  by  a 
one-dimensional  turbulent  model.6  It  may  be  noted  that  solutions  for 
this  set  of  graphs  were  obtained  by  using  the  main  stream  wind 
velocity  U  equal  to  10  m/s.  To  investigate  effects  of  cloud-top- 
warm-air  entrainment  on  stability  of  the  cloud  layer  shown  in  figure 
2,  a  stream  function  shown  in  figure  3  will  be  superimposed  to  the 
main  flow.  In  generating  the  stream  function  shown  in  figure  3,  the 
top-down  entrainment  was  assumed  to  be  at  a  maximum  rate  of  6 
cm/s  is  at  the  top-left  comer  and  the  rate  was  reduced  sinusoidally 
to  zero  at  the  top-right  comer.  In  addition,  it  is  assumed  that 
potential  temperature  of  the  entrained  air  is  at  five  degrees 
centigrade  higher  than  that  of  the  cloud  at  the  top. 

Using  the  initial  conditions  and  entrainment  rates  described 
above,  the  present  large-eddy  simulation  computer  program  has 
been  run.  Shown  in  figure  4  are  contours  of  the  initial  steady-state 
potential  temperature  values,  total  moisture  mixing  ratio  values  and 
liquid  moisture  Mixing  values,  respectively.  Predicted  dynamic 
responses  of  the  cloud  layer’s  liquid  moisture  content  to  the  warm- 
air  entrainment  at  the  top  are  shown  in  figure  5.  From  snapshots 
shown  in  this  figure  of  the  cloud  contours  at  different  times  (i.e.,  at 
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eight  Hemet  eight  ilometers  eighfeitomet  eighfeitomet 


@  Half  hour  later 


Horizontal  distance  in  kilometers 


Fig.  5:  Predicted  Effects  of  Warm-Air  Enrtrainment  on  the  Cloud  Stability 
(Plotted  contours  at  different  time  are  for  liquid  water  humidity  ratio 
values  in  grams  per  kilogram) 
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